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Progress Report for the Period July 1991 to December 1991 
NASA LaRC Grant No. NAG-1-1215: 

Improvements to a Method lor the Geometrically Nonlinear Analysis ol 
Compressively Loaded Stittened Composite Panels 



This report describes progress made during the period July 1991 to December 1991 on 

the tasks identified in the technical proposals for the subject grant. The plans lor further effort 

on each of the tasks are outlined. The computer implementation of the method of analysis 

under development is referred to In this document as NLPAN. o 
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i) Implementation of Conti nuation Methods <> 


in the progress report for Ihe preceding period of performance Ml. an outline was given 
of the theory behind the implementation of two so-called "continuation methods, namely Rlks 
method 121 and Thurston's method [31, within the NLPAN computer code. The methods are 
used for traversing limit points and bifurcation points, respectively, in the equilibrium 
load/response behavior of stiffened panels. Problems were encountered in the present re- 
porting period which led to some changes both to the implementation of the two continuation 
methods, and to the fundamental formulation used in NLPAN. as described In Ref. [41. These 
problems and changes are brier, described here. Complete details of the modified theory 
are included in a theory document which is under preparation, and which will be releaset 

during the next period of performance. 

The set of nonlinear algebraic equations which are ultimately solved in order to compu, 
equilibrium solutions are written symbolically as 
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where lambda is the generalized load parameter, and qr is the r* modal amplitude 
, , , , 2 M. The following two matrices play roles in the two continuation methods 
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and in Ref [1] if was stated that the tangent stitfness matrix. [00. is symmetric, and that the 
negative matrix [-0*] is both symmetric and positive dehnite. Both of these assumptions 
have been found to be false. If a function f, is derived from a total potential energy function, 

7 T, by the relationship 


f i = 


dn 

dq. 


(4) 


then the corresponding tangent stiffness defined by egua.ion (2) will indeed be positive defi- 
nite [2] it was thus discovered that the virtual work approach behind NLPAN as described in 
,4, does no,, as was prevtousi, believed, correspond ,0 a statement of stationary tota, poten- 

tial energy. 

Because of the desire to have a symmetric tangent stitfness matrix for the sake of the 
solution procedures for the nonlinear algebraic equations, the direct application o, the princt- 
p,e o. virtual work was abandoned and replaced with an approach which begins w„h a total 
potentia, energy expression The new approach has led to certain restrictions regardtng the 
selection ot load control versus displacement control tor the ln-p,ane loads on a pane,, and 
these restrictions are summarized in Table 1. 

The assumption tha, the matrix [-0*] is symmetric and positive-dehnite was tunda- 
men, a, to the approach descrtbed in Ret. ,1, tor monitoring the stability o, the equilibrium 
state, and because the assumption was incorrect, a change was required. The g 

problem posed in [1], 


has 


([D°] + ^[D a ])(^} = { 0} . *“1.2 M 

been replaced with by the following eigenvalue problem: 
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where [/] is the identity matrix. This eigenvalue problem is used to determine the stability of 
an equilibrium solution; when the equilibrium path approaches a critical stability state, the 
lowest eigenvalue, a>u goes to zero. 


Results computed using both Riks' method and Thurston's method are presented in Fig- 
ure 1. The configuration modelled is a rectangular isotropic plate with bi-axial loading applied 
through control of the in-plane normal displacements at each edge. Other aspects of the plate 
configuration are given in the figure. Two-mode analyses were performed for the purpose of 
illustrating the operation of the two methods discussed here. First, a perfect plate was mod- 
elled, and Thurston's method was used for navigating past the three bifurcation points, labeled 
A, B, and C, which are encountered along the equilibrium path, labeled Path 1. Next, an 
imperfect plate was modeled, where the imperfection shape has contributions in the shapes 
of both of the mode shapes used. The equilibrium path of the imperfect plate, labeled Path 
2, exhibits two limit points, labeled D and E, and these limit points were traversed in the 
analysis using Riks' method. Both of the curves were generated by the same automated 
computer procedure which integrates the two continuation methods. There are two dots on 
each equilibrium path. One dot on each path indicates where the equilibrium state of the plate 
becomes unstable for a load steadily increased from zero, and the other dot on each path in- 
dicates the state at which an actual plate would come to rest after a dynamic snap. 

One difficulty with the two solution methods discussed here is that for some configura- 
tions, the second eigenvalue of equation (6) goes to zero, indicating, it is believed, that there 
are multiple unstable equilibrium paths which intersect. No robust logic for tracking through 
such a tangle has been established, and NLPAN has been programmed to terminate the sol- 
ution process when this situation is encountered. The dynamic analysis discussed below is 
designed to circumvent this complication. 

jit Dynamic Analysis Capability 

The purpose of the dynamic analysis capability is to analytically locate the stable equi- 
librium state sought by a structure which has passed from a stable equilibrium state to an 
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unstable one, and which thus begins a transient dynamic response. It is assumed that the 
dynamic response occurs at a fixed value of the generalized load, namely the value at which 
the equilibrium solution becomes unstable, and that damping is present to dissipate energy. 
The detailed formulation of the dynamic analysis is provided in a theory document which will 
be released during the next reporting period. A brief summary of the approach is given here. 

The form of the nonlinear algebraic equations governing equilibrium is given by 
(AC, L + C°) + q/ACjj 4- C°) + + C° k ) + q^q^AC^ + C°^) = 0 , / = 1. 2. ... ,/M (7) 

where summation is implied over repeated indices, A is the generalized load parameter, and 

the M modal amplitudes, q, (/= 1,2 M), are the generalized coordinates of the structure. 

The above equations represent the condition of stationary total potential energy. Hamilton s 
principle is applied in order to obtain the equations of motion, and an equation of the following 
form is ultimately obtained: 

qfjj + qpij + QjC/j + qflkCfjk + qflifls + C l]kf = 0 , / = 1 , 2 M (8) 

where q y and q y are the first and second time derivatives, respectively, of the modal ampli- 
tude, and where 

C, = A* C ( L + Cf C^TC^+c; (g) 

Cjj = /l* cjj + C,j Cjjkt — A CjjM + C,j k f 

where A* is the value of the generalized load parameter at the critical equilibrium state, and 
P / 

d <j = E 

p = i V 4 


(10a) 



(10b) 
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where p is a plate strip index number, P is the total number of plate strips in the structure, p 
is a viscous damping coefficient with units of force per unit area, per unit velocity, m is the 
mass per unit area of a plate strip, and v, and w, are components of the buckling modes. 

The Newmark direct integration procedure [5] will be used to solve equations (8) at a 
series of uniform time steps. Because the goal of the analysis is to locate a new equilibrium 
solution, the selection of the viscous damping coefficient and the time increment need not be 
made with the concern of obtaining exact or realistic dynamical solutions. Nonetheless, the 
selection of the two parameters mentioned will effect the performance of the dynamic solution 
strategy. The basis for the selection of the parameters has not yet been established, but the 
issue will be taken up when coding of the dynamic analysis procedure is complete. 

iii) Additional Boundary Condition Options for the Panel Ends 

The formulation has been completed for two new varieties of boundary conditions which 
can be imposed at the longitudinal ends of a panel, namely clamped ends, and eccentrically 
applied end loading. For the clamped-end case, the condition which is sought is that the axial 
displacement is the same for all points at each end of the panel. Because this condition can 
not generally be satisfied exactly, the error in the satisfaction of the boundary condition is 
minimized in the least squares sense. This is accomplished through the use of the penalty 
function method (6], 

In modelling load eccentricity, the goal is to simulate a line across the end section of the 
panel along which the load is assumed to be applied in the fashion of a knife-edge. The 
modelling of this support condition is more complicated than for the clamped-end simulation. 
Several points in the panel cross-section are specified which define the line of load applica- 
tion. These points are constrained to move together using the penalty parameter method, and 
an additional variable parameter is added to the expression of displacements. The new pa- 
rameter controls an end-shortening correction to the first-order displacements, in order to 
provide the condition that the pivoting (due to the first-order displacements) of the panel ends 
about the eccentric load lines occurs without any effective change in the end-shortening value. 
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Coding of the method for modelling the two different end-support conditions mentioned 
above in in progress, and the details of the formulation are in a theory document to be re- 
leased in the next reporting period. An additional task planned is to develop and implement 
in NLPAN the ability to model elastic restraint at the panel ends. 

iv) Transverse Pressure Loading 

As described in the previous progress report, the ability to model transverse pressure 
loading has been developed and encoded in the NLPAN computer program. Further progress 
during the present reporting period has been the implementation of a method for controlling 
asynchronous application of the in-plane loading and the transverse pressure loading. De- 
noting the in-plane load parameter as A and the pressure load parameter as /?, a series of 
target load pairs, (A„, /?,) = (0,0), (A,, /?,), (A 2 , /? 2 ), ... is specified which bound a series of load 
ranges. In each load range, parameters X and vary linearly with a new, generalized load 
parameter, (. where £ varies from zero to one within each load interval. Thus, the two load 
parameters are replaced with a single one, and the solution strategies for the nonlinear al- 
gebraic equations can be readily applied. It is noted that the in-plane load parameter, A, may 
control either in-plane normal loads or in-plane normal displacements, consistent with the 
selections made from Table 1. The details of the theory for transverse pressure loading and 
the method of controlling asynchronous loading are contained in a theory document to be re- 
leased in the next period of performance. 

As an example of the application of asynchronous loading, NLPAN was used to model a 
long isotropic plate subjected to combined axial compression and transverse pressure. The 
configuration is one for which results were computed by Levy et al [7], and details of the 
configuration are provided in Figure 2a). For verification purposes, four buckling mode shapes 
were incorporated in the NLPAN analysis corresponding to the four double sine functions used 
in the Ref [7] to represent out-of-plane displacements. The computed values of axial load 
versus end-shortening are presented in Figure 2b), and it can be seen that there is essentially 
exact agreement between the NLPAN predictions and the results of Ref. [7], The NLPAN 
analysis was conducted using two load intervals, the first of which corresponded to the ap- 
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plication of transverse pressure with zero net end load. It can be seen that the panel shortens 
in length due to the application of the transverse pressure. The second interval simulates the 
application of axial loading while the transverse pressure maintained at a constant value. 
Rik's method was called on by the NLPAN analysis in order to traverse the multiple limit 
points which are encountered. During the next period of performance, NLPAN will be used to 
model stiffened panel configurations subjected to pressure loading. 

v) Second-Order Displacement Fields 

As discussed in the previous progress report, there is a discrepancy in the functional 
form of the second-order displacement fields with regard to satisfaction of the boundary con- 
ditions at the longitudinal ends of a panel. The transverse displacement components of the 
buckling modes, v, and w„ vary as sin(mrrx/L) in the longitudinal direction, where m is the in- 
teger halfwave number, and thus these displacement components go to zero at the longitudi- 
nal ends of a panel. In contrast, the transverse displacements of the second-order 
displacement fields, v H and w, h vary as cos (mn-x/L) in the longitudinal direction, where m is 
an integer, and thus w, and w, ; do not in general go to zero at the panel ends. 

Recent efforts toward addressing this problem have focused on the fact that the ordinary 
differential equations governing the y-dependence of the second-order displacement functions 
are identical to the ordinary differential equations which govern the y-dependence of the 
buckling eigenfunctions, with the distinction that the latter differential equations are homoge- 
neous, whereas the former equations are not, and where it is understood that the two sets of 
ordinary differential equations referred to depend on the the associated halfwave number, m 
or m, respectively. The similarity of the two sets of differential equations has several impli- 
cations. First, the second-order displacement fields contain contributions which correspond 
to buckling mode shapes, except that the phase difference with respect to the longitudinal di- 
rection causes the points of maximum displacement to occur at the ends of the panel, where 
they violate the boundary conditions. Second, the degree to which particular buckling mode 
shapes are present in a second-order displacement field depends on the reference value of 
the load parameter used in computing the second-order field. The shape of the field is thus 
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highly dependent of the reference load level selected, and amplitude of the field becomes 
unbounded when the reference load level coincides with a buckling eigenvalue for the partic- 
ular halfwave number associated with the Held being computed. 

It would seem that the ideally, the second-order fields should be forced to be orthogonal 
to all buckling mode shapes, in which case the second-order displacement fields would truly 
be second-order in nature in that they would not duplicate displacement contributions which 
could be represented by the inclusion of buckling mode shapes in the analysis. It is suspected 
that the second-order fields would then be both load-independent, and would be much more 
in accordance with the boundary conditions at the panel ends. In analyses of thin-walled 
structural sections, Benito and Sridharan [8] enforced orthogonality of the second-order dis- 
placement fields with respect to the buckling modes through the use of the Lagrange multi- 
plier method. This approach has the disadvantage of eliminating the banded quality of the 
coefficient matrix of the linear system of equations which is solved to obtain a discrete rep- 
resentation of a second-order field by the current method. Thus, enforcing the orthogonality 
condition by this method may add significantly to the computational cost. This and other 
means of enforcing the will be investigated further during the next period of performance. 

vi) Additional Tasks 

Additional efforts are planned in the next reporting period, consistent with the continua- 
tion proposal recently submitted to NASA Langley Research Center. The most significant of 
these tasks is to develop the ability within NLPAN to model the effects of thermal loading. 
An additional task is to improve the computational economy of NLPAN, both with respect to 
central processing unit execution times, and with respect to computer core memory require- 
ments. Reductions in the core memory requirements will result in an increase in the modelling 
accuracy available, because of the ability to incorporate a larger number of mode shapes in 
an analysis. Finally, an attempt will be made to automate the process of selecting the mode 
shapes to be incorporated in an NLPAN analysis. Currently, this selection is left to the dis- 
cretion of the program user, but systematic approaches have been developed for selecting 
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mode sets for certain types of configurations, and it is planned to implement these approaches 
in an automated fashion. 

vii) Results for an 1-Stiffened Panei with a C omplex Cross-Section 


NLPAN was used to model the postbuckled response of a stiffened graphite/epoxy panei 
with a complicated cross section, which called for the use of offsets between the plate strip 
edges and their associated node lines. The ability of NLPAN to account for these offsets has 
not previously been exercised. The configuration modelled is a stiffened panel for which ex- 
perimental results are given by Starnes et al [9], The panel has four stiffeners, and was flat- 
end tested in uniaxial loading. Some details of the panel configuration are given in Figure 
3a-b), and complete details are given in Ref. [9], A unit cell representation was used in the 
NLPAN analysis, and a schematic drawing of the cross-section of the unit cell is shown in 
Figure 3c). The unit cell model employed 17 plate strips and 12 node lines, and a three-step 
representation of the tapered flanges of the test article was used. The node-line offsets re- 
ferred to above are visible in Figure 3c. Symmetry conditions were enforced at the side-edges 
of the unit cell. A four-mode NLPAN analysis was used, where the dominant mode was the 
primary local-buckling mode having five halfwaves along the length of the panel, and the three 
additional mode shapes were selected based on their ability to provide the refinement of the 
displacement field (with respect to the primary buckling mode shape) which takes place in 
postbuckling. An initial imperfection was used in the NLPAN model in the shape of the pri- 
mary buckling mode and having an amplitude of one percent of the panel skin thickness. 


Results for end load versus end shortening are presented in Figure 4a), and the agree- 
ment between the NLPAN predictions and the experimental results is fairly good, though 
NLPAN predicts slight'y less axial stiffness in postbuckling than was observed experimentally. 
(The NLPAN model used the unsupported length of the panel for the postbuckling analysis, 
with a correction added to account for the compression of the panel ends within the potted 
ends.) The distribution of longitudinal membrane strains (in the panel skin) across the width 
of the center bay at the mid-length of the panel are shown for three different load levels in 
Figure 4b). The agreement between the theory and the experiment is good except for a slight 
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over-prediction of strain levels by NLPAN near the center of the panel at the higher load lev- 
els. 

In Figure 5, results are presented for longitudinal strains on opposing surfaces at two lo- 
cations on the mid-length of the panel, as affected by the panel loading. Figure 5a) shows the 
location of strain gages used to obtain the experimental results. The strain levels at the center 
of the panel skin are plotted in Figure 5b), and the agreement between the experiment and 
theoretical results is fairly good. NLPAN over-predicts the strain amplitudes slightly at the 
higher load levels, and this is consistent with the discrepancy noted regarding the results of 
Figure 4b). The strain levels at a stiffener flange location are plotted in Figure 5c). The strains 
predicted by NLPAN run uniformly about 10% less than the reported experimental values. 
There was some disagreement between the slopes of the curves in prebuckling, so the ex- 
perimental results were plotted a second time with the strain values scaled down so that the 
prebuckling slopes of the theory and experiment match; the modified experimental curves lie 
nearly on top of the theoretically curves. The reason for the discrepancy in the prebuckling 
slopes is not apparent, but the qualitative agreement between the theoretical results and the 
experimental results provides some confidence that the node-line offset feature of the NLPAN 

program is working properly. 
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Control Param.: 
L- Edge Load 
D- Edge Displ. 
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Component 2 
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Imperfect plate: 







a/b=4 

Simple support on all edges 



Assumed form for w(x,y): 
w(x,y)=[w 1 1 sin(l n x/a ) +w 31 sin(3 rc x/a ) 

+wg-| sin(5 7t x/a ) +wy-j sin(7 it x/a ) ] sin( 7t y/b ) 



Figure 2. Isotropic rectangular plate under combined axial and pressure loading. 
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b) Membrane strains in the skin, mid-length, center bay. 


igure 4. Analytical and experimental results for a stiffened graphite/epoxy panel in postbucklmg 
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a) Strain gage placement at the mid-length of the center bay. 



b) Strains at the center of the skin. 



c) Strains at the stiffener flange. 


Figure 5. Longitudinal surface strains on a stiffened graphite/epoxy panel in postbuckling. 
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